Dynamic critical behavior of the worm algorithm for the Ising model 
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We study the dynamic critical behavior of the worm algorithm for the two- and three-dimensional 
Ising models, by Monte Carlo simulation. The autocorrelation functions exhibit an unusual three- 
time-scale behavior. As a practical matter, the worm algorithm is slightly more efficient than 
Swendsen-Wang for simulating the two-point function of the three-dimensional Ising model. 
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Monte Carlo simulations in statistical mechanics [l| 
and quantum field theory [2J typically suffer from criti- 
cal slowing- down 0, 0| : the autocorrelation (relaxation) 
time r diverges as the critical point is approached, most 
often like r ~ £ z , where £ is the spatial correlation length 
and z is a dynamic critical exponent. For conventional lo- 
cal algorithms, one usually has z 2. This effect severely 
limits the efficiency of Monte Carlo studies of critical phe- 
nomena in statistical mechanics and of the continuum 
limit in quantum field theory. 

One approach to circumventing this slowing-down in- 
volves replacing the underlying spins or fields by an al- 
ternate representation, obtained from the original model 
by algebraic transformation. For instance, the celebrated 
Swendsen-Wang (SW) cluster algorithm [5( simulates the 
q-state ferromagnetic Potts model 0, 0] by passing back 
and forth between the Potts spin representation and the 
Fortuin-Kasteleyn bond representation [1, 0|. In this 
Letter we shall stud y a nother such algorithm, namely 
the worm algorithm |lOl . [H} , which simulates the high- 
temperature graphs of the spin model, considered as a 
statistical- mechanical model in their own right. Surpris- 
ingly, no systematic study of the dynamic critical behav- 
ior of the worm algorithm has heretofore been carried 
out, even in the simplest case of the Ising model. As we 
shall show, the worm algorithm presents some unusual 
dynamical features, which combine to make it an extraor- 
dinarily efficient algorithm for simulating some (but not 
all) aspects of the three-dimensional Ising model. Indeed, 
it is surprising (at least to us) that an algorithm based 
on local ( "worm diffusion" ) moves could perform so well. 

We consider the zero-field ferromagnetic Ising model, 
with nearest-neighbor coupling J > 0, on a connected fi- 
nite graph G = (V, E) with vertex set V and edge set E. 
The high-temperature graphs of this model are subsets 
A C E of "occupied bonds" . For any bond configuration 
A C E, we denote by dA the set of vertices that touch 
an odd number of bonds of A. We focus attention on the 
set 6>0 = {A: dA = 0} of "vacuum graphs" and the sets 
S x ,y = {A: dA = {x,y}} of "two-point-function graphs", 
with the convention that <S X!X = 6>0. Clearly, Sp.,, = S y ^ x . 
The standard high-temperature expansion 12] states 



that Z = E A es ™ lAl and ZG(x,y) = T,Aes x>y wlA{ > 
where Z is the Ising partition function (up to an uninter- 
esting prefactor), G(x,y) = (o~ x a y } is the Ising two-point 
correlation function, w = tanh J, and \A\ is the number 
of occupied bonds in the configuration A. 

Instead of simulating the space of Ising spin configu- 
rations, the worm algorithm simulates a space of high- 
temperature graphs. More specifically, the configuration 
space S of our version of the worm algorithm consists of 
ordered triplets (A, x, y) with x, y € V and A 6 S x>y , i.e., 
A is a bond configuration having odd degree at x and y 
(unless x — y) and even degree at all other sites. We set 
the weight of a configuration (A, x, y) to be F(x,y) 
where F is an arbitrary nonnegative function. It follows 
that the "partition function" of our ensemble is 



Z = 



E ' 

{A,x,y)£S 



l F(x,y) = Z 



E 

x,y£V 



F(x,y)G(x,y) 



(1) 

We shall usually take F = 1, so that Z = Z(M 2 ), where 
M. = a x is the total magnetization; in a translation- 
invariant situation Z = ZV\, where V is the volume and 
X is the Ising susceptibility. This configuration space is 
clearly tailored to studying the Ising two-point function. 

The elementary move of the worm algorithm is as fol- 
lows: Pick uniformly at random one of the two "end- 
points" (say, x) and one of the edges emanating from x 
(say, e = xx'). Propose to move from the current con- 
figuration (A, x, y) to the new configuration (AAe, x' , y), 
where A denotes symmetric difference (i.e., delete the 
bond e if it is present, or insert it if it is absent). Then ac- 
cept or reject this move according to either the Metropo- 
lis or the heat-bath criterion. For instance, in the heat- 
bath version, the configuration with (resp. without) e 
gets probability w/ (1 + w) [resp. 1/(1 + w)]. This simu- 
lates the distribution |T]) with F(x,y) = d x d y , where d x 
is the degree of the vertex x in G (i.e., the number of 
nearest neighbors). Other choices of F can be simulated 
by an appropriate Metropolis accept-reject step. 

Additional moves can optionally be added. Because of 
the symmetry x <-> y, we can interchange x and y with 
probability 1/2 after each worm move. More interest- 
ingly, whenever we reach x = y, we can move the end- 
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points from (x, x) to a randomly chosen (x 1 , x 1 ). Finally, 
we can add "local" moves A i— > AAB (to be accepted or 
rejected according to the Metropolis criterion), where B 
is any bond configuration having DB = 0, e.g. a plaque- 
tte or a winding cycle. 

We remark that the "worm" idea is very general: en- 
large a state space of "vacuum" (Eulerian) bond configu- 
rations to include a pair of "dislocations" , and then move 
those dislocations by random walk. This idea can be ap- 
plied, in particular, to the hexagonal-lattice 0{n) loop 
model [13| at general n, and to vertex models. 

In this Letter we report detailed measurements of the 
dynamic critical behavior of the worm algorithm and 
some of its variants, for two- and three-dimensional Ising 
models at criticality, on lattices of size L d with periodic 
boundary conditions. We shall measure time in units of 
"hits" of a single bond, but we stress that the natural 
unit of time is one "sweep" of the lattice, consisting of 
L d hits. We write z = x — y for the vector distance 
between the endpoints, and we define the observables 



and Tr, 



e zpz . In our simulations we mea- 



sured the number AT = \A\ of occupied bonds as well as 
the short-distance observable T>q and the low-momentum 
observable T\ ovl — (l/2d)J2\ p \ = 2w/L-^p- Please note that 
(Do) = l/x, {Km) = G(p)/x for |p| = 2n/L, and (TV) = 
w d/dwlog(Zx)- In particular, the second-moment cor- 
relation length H is C = ((^low) -1 - l) 1/2 /[2sin(vr/ J L)]. 

For any observable O, let po(t) be its normalized auto- 
correlation function in the stationary stochastic process. 
Then define the exponential autocorrelation time 



Texp.o = limsup- 



1*1 



t^±oo - log |/30 (t) I 

and the integrated autocorrelation time 



(2) 



(3) 



t=-c 



Typically all observables O (except those that, for sym- 
metry reasons, are "orthogonal" to the slowest mode) 
have the same value t CXPi o = r cxp . However, they may 
have very different amplitudes of "overlap" with this 
slowest mode; in particular, they may have very differ- 
ent values of the integrated autocorrelation time, which 
controls the efficiency of Monte Carlo simulations [4(. 
We define dynamic critical exponents z cxp and z- mtt o by 
T ox P ~ £ Zcxp and Tint,© ~ £ Zint '°, where time is measured 
in "sweeps" . On a finite lattice at criticality, £ can here 
be replaced by L. 

Before presenting our numerical results, let us make 
some heuristic predictions for the dynamic behavior of 
the worm algorithm. Suppose first that the bond config- 
uration A is at all times completely equilibrated for the 
given endpoints x,y [l5| . Then z = x — y performs a 
random walk with drift having equilibrium distribution 



G(z)/x- In the simplest case G = 1 (which corresponds 
to the zero-temperature limit J = oo), the eigenvectors 
of this random walk are J- p , with eigenvalues (in the 
heat-bath version) A p = [1 + (1/d) Yli=i cospJ/2. In 



' cxp 



L . Furthermore, the autocorrelation 



particular, 

function of 2? a is const x X) P7 £o ■ ^ n ^ ne limit L — > oo, 
this tends to -Pt(O), the return probability of the corre- 
sponding random walk on Z d ; in particular, for large t 
it behaves like t~ d / 2 . It follows that Tj nti x> a ~ L 2 ~ d in 
d < 2, logL in d = 2, and L° in d > 2. We expect 
these scaling predictions to hold near the critical point 
in d = 1 (since J c = oo) and in the low-temperature 
regime (J > J c ) in d > 2. A more complicated behavior 
is likely to occur in the critical regime in d > 2, where 



G(z) 



-(d-2+7j) 



as |z| — > oo. A Fokker-Planck analy- 



sis under the hypothesis of perfect equilibration of bonds 
predicts p Pa (t) ~ H 1 " 1 ^ 2 ). 

It should be noted, however, that Vq estimates x vla - a 
"rare" event, i.e. a binomial random variable with proba- 
bility 1/x- The variance of this random variable is also of 
order l/x, so that ~ x samples are needed to get a rela- 
tive variance of order 1. This is an example of "statistical 
inefficiency due to large static variance" . 

Note, finally, that a simple variational argument (fol- 
lowing the model in 0]) shows that Zint.jv" ^ ex./ v. 

We began by simulating the worm algorithm at the 
critical point in d = 1 and in the low-temperature phase 
in d = 2,3. The autocorrelation functions of J-\ ow and Vq 
behave exactly as predicted. The autocorrelation func- 
tion of M is essentially a pure exponential, with autocor- 
relation time ~ L 2 in d = 1, ~ L 2 \ogL in d = 2, and 
- L 3 in d = 3. 

We next simulated the square-lattice (d — 2) Ising 
model at criticality (w c = \[2 — 1) on lattices 4 < L < 
2048. The simulation lengths varied from 5 x 10 11 hits 
(L = 4) to 5.4x 10 14 hits [L = 2048). These runs used ap- 
proximately 6.9 yr CPU time on a 3.2 GHz Xeon EM64T 
processor. Estimates of the static quantities Xi £ an d 
E agreed within error bars with previous high-precision 
simulations using the SW algorithm [lij ]. 

The slowest mode is the number M of occupied bonds, 
and the decay of pj\f{t) is very close to pure exponen- 
tial (Fig. HJ. A fit Ti n t,AA/L 2 - AL Z (resp. AL Z + B) 
yields z cxp = z- lnt ,M ~ 0.379 (resp. 0.338). But a behav- 
ior T- lnt _^f /L 2 ~ log 2 L is also conceivable. 

A more complicated behavior is exhibited by ^low: 
contrary to our simple-minded prediction, its decay is 
far from a pure exponential. Rather, it appears to scale 
like pjc low (i) w f(t/L d+z ) with a scaling function / 
that shows an initial power-law decay f(x) ~ x~ r with 
r ps 3.05 but then bends towards an unknown smaller 
power (Fig. O. We find z' w 0.315. It is not clear 
whether z' equals z cxp or is slightly smaller. 

The most interesting behavior of all is shown by T>o, 
which exhibits significant decorrelation on a time scale of 
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Figure 1: Autocorrelation function pAf(i) versus t/n n tM for 
the critical two-dimensional Ising model. 
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Figure 2: Autocorrelation function pr la „(t) versus 1 + 
at/L d+z ' with a = 2.5, z = 0.315 for the critical two- 
dimensional Ising model. Black line shows initial decay with 
power r ~ 3.05. 
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Figure 3: Autocorrelation function pT> (t) versus t for the 
critical two-dimensional Ising model. Asymptotic straight line 
has power s — 0.75. Dashed line has power 1 — rj/2 = 7/8. 
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Figure 4: t s pT> (t) versus i/rint.A/ - f° r the critical two- 
dimensional Ising model, with s — 0.75. Dashed line shows 
the decay rate r eX p. 



order 1 hits (Fig. [3]). The data fall in the limit L — > oo on 
a beautiful scaling curve PT> (t) — g{t), with g{x) ~ x~ s 
as x — > oo. We find s w 0.75, which is clearly smaller 
than our prediction s = 1 — r\j1 = 7/8 based on perfect 
equilibration of bonds. Apparently this latter prediction 
provides a lower bound on pv (t) an d hence an upper 
bound on the exponent s. 

Finally, we analyzed the universal crossover from 
short-time to long-time behavior in px> , which we hy- 
pothesize is of the form pv Q {t) — g(t) h(t/L d+Zoxp ), by 
plotting t s p-D {t) versus i/ T int,7V- A fairly clear scaling 
curve is seen (Fig. [4]), though it is noisy for large lattices. 
Using this scaling Ansatz to compute the area under the 
curve of px> (t), we conclude that 

_ / -sd + (1 - s)z cxp if s < 1 f .s 
Zint < V ° \-d if S >1 (4) 

Using z cxp m 0.338 and s sa 0.75, we find z- m t.v ~ —1-42. 



We also studied the variant of worm algorithm with the 
move (x, x) — > (x' , x'). We find that, in the limit L — > oo, 
this move has no effect (i.e., the ratio of autocorrelation 
times with and without the move tends to 1). Evidently 
the diffusion of the endpoints x, y around the lattice in 
the basic worm algorithm is already sufficient. 

We next simulated the simple-cubic (d — 3) Ising 
model at the estimated critical point J c = 0.22165455 
on lattices 4 < L < 256. The simulation lengths var- 
ied from 5 x 10 10 hits [L = 4) to 7.2 x 10 14 hits [L = 256). 
These runs used approximately 6.6 yr CPU time. 

Once again the slowest mode is A/", and the decay is 
very close to pure exponential. In contrast to d = 2, 
the exponent z cxp = Zmt.A/' now appears to equal a/v w 
0.174. Again T\ ow exhibits a scaling curve with a very 
strong bending and an unknown asymptotic decay ex- 
ponent r, with z' w —0.15 < z cxp (Fig. [5]). Finally, 
T>q exhibits a clear scaling with exponent s w 0.66 
(Fig. [5]). Again this exponent is smaller than our pre- 
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Figure 5: Autocorrelation function pr lo „{i) versus 1 + 
at/L d+z with a = 1.0, z = -0.15 for the critical three- 
dimensional Ising model. 



d = 2 but significantly better in d = 3, where zgw ~ 0.46 
pjj . For estimating £ using .Flow, the exponent Zi n t.r lom is 
uncertain because of the uncertainty on the decay expo- 
nent r (Figs.[2]and[5|), but it must lie somewhere between 



and 



i.e. from 0.31 to 0.38 in d = 2 and from —0.15 



to 0.17 in d = 3. This is again slightly worse than SW in 
d = 2 but better in d — 3. We conclude that the worm 
algorithm is, asymptotically as I — > oo, the most effi- 
cient algorithm currently available for simulating x an d 
£ in the three-dimensional Ising model. In practice, our 
data show (l6j that the worm algorithm outperforms SW 
when L > 32, at a rate that grows like L~ 32 . 

Details of these simulations and their data analysis will 
be reported separately [ijj]- 
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Figure 6: Autocorrelation function pT> (t) versus t for the 
critical three-dimensional Ising model. Asymptotic straight 
line has power s = 0.66. Dashed straight line has power 
1 — 77/2 = 0.982. 

diction s = 1 — 77/2 w 0.982 [13] based on perfect equili- 
bration of bonds. From ((4]) we obtain z- m t.T> ~ —1.92. 

In summary, the worm dynamics for the critical Ising 
model appears to exhibit decorrelation on three different 
time scales: Af has an almost-pure-exponential decay on 
the very long time scale L d+Z " x p (in "hits"); JF low has 
a complicated power-law decay on the long time scale 
L d+Z ; and T> has a simple power-law decay on the short 
time scale L°. This is analogous to but more complex 
than the two-time-scale behavior recently observed in the 
Sweeny dynamics for the random-cluster model [ljj ]. 

The practical efficiency of the worm algorithm de- 
pends on the observable. For estimating x using Vq, 
the autocorrelation time Ti nt .v ~ L Zint l3 o [cf. (g])] must 
be multiplied by a factor x ~ due to static vari- 

ance, leading to an "effective dynamic critical exponent" 
ZeS,v = <Zint,x> + ~ 0.33 in d = 2 and « 0.04 in 

d = 3. This is slightly worse than the SW algorithm in 
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